Computational Consequences of Neglected First-Order van der Waals Forces 
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We have computed the widely neglected first-order interaction between neutral atoms. At inter- 
atomic separations typical of condensed media, it is nearly equal to the — f/r 6 second-order London 
energy inferred from interactions in gases. Our results, without the exchange forces that lead to cova- 
lent bonding, suggest that the quality of non-bonding attraction between neutral atoms of molecules 
in condensed media differs from the — 1/r 6 form usually ascribed to it. If we add first-order and all 
second-order terms, we obtain energies nearly five times the —1/r 6 London energies which dominate 
only at the atomic separations found in gases. For computation, we propose a practical, accurate 
form of energy to replace the qualitatively inaccurate Lennard- Jones and harmonic forms casually 
assumed to hold at the interatomic separations found in condensed media. 

PACS numbers: 82.35.Pq,34.20.Cf,77.84. Jd 



Does the van der Waals potential energy between two 
neutral atoms fall off as — 1/r 6 at separations of a few A? 
Is the widely neglected first-order van der Waals potential 
zero at the interaction distances within condensed media? 
Should one use a Lennard- Jones 6-12 potential in com- 
puter simulations? Because it is no longer necessary to 
sacrifice accuracy for computational convenience, we nu- 
merically re-examine these old questions, find that com- 
mon practice leads to qualitative error, and propose an 
accurate but practical alternative to the mathematically 
convenient but potentially misleading Lennard- Jones (L- 
J) and harmonic forms. 

Motivated by an intriguing footnote in Landau and 
Lifshitz we compute the first-order van der Waals 
potential for two hydrogen atoms, without the exchange 
interaction that creates covalent bonds. It is nearly equal 
to the 1/r 6 second-order London energy found from the 
same wave functions when evaluated over a range of sep- 
arations where both exceed kT at room temperature. We 
then compute the potential energy to second-order, ob- 
taining an expression that reassuringly reproduces the 
expected dominance of the 1/r 6 London term as r — > oo. 
The net result at separations 3 ciq < r < 5 ao (Bohr ra- 
dius ao) is an energy 4 to 5 times that expected from a 
London interaction, which emerges to dominate only in 
the r > 10 ao limit. 

We next examine the additivity of the interactions of 
several particles, usually taken to be additive in first and 
second order; the Axilrod-Teller-Muto force @, Q is a 
third-order effect. In fact, energies add in first-order per- 
turbation theory but not in second order unless one ne- 
glects first-order terms. 

Constructed to fit the measured equilibrium separation 
r s and the measured condensation energy energy V(r s ) 
as well as to have the large-r asymptotic form 1/r 6 , the 



Lennard-Jones potential is too hard for r < r s and varies 
incorrectly for r > r s . Harmonic potentials work only 
near energy minima. We propose a simple, phenomeno- 
logical potential that can fit spectroscopic data and that 
behaves appropriately at all distances. 

For our calculations, we take the hamiltonian for two 
interacting hydrogen atoms to be H = Hq + W in which 
Ho is a sum of two isolatcd-hydrogcn-atom hamiltonians 
and the perturbation IF is a function of the distance r 
between the two protons 
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Here r = \r\ and e is the charge of the electron in units 
with a = e 2 /(hc) w 1/137. The first-order correction 
to the energy is AE\(r) — (0|IF|0) in which the unper- 
turbed ground state |0) = 1 100, 100) is two Is states sep- 
arated by a distance r. Its wave function is (fi)7^|0) = 
(l/TOp) exp(— ri/oo — 7-2/00) where ao = h 2 /(me 2 ) is 
the Bohr radius, r\ = ri|, and r 2 = t^I- Integration of 
(0\W\0) yields 
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London and others 0, @, @1 applied second-order pertur- 
bation theory to the potential ([1]), expanded for large r, 
and got AE L (r) = -6.50 mc 2 a 2 a^/r 6 . 

Fig.Q]shows that the two energies AEi(r) and AE L (r) 
are nearly equal for 3ao < r < 5a , where both are reli- 
able and greater than kT = 1/40 eV. Because AEx{r) is 
the first-order result, and AEr.(r) is an approximation to 
second-order perturbation theory, the two can be added. 
That is, to (approximate) second order in the perturba- 
tion W and within the London approximation, the van 
der Waals energy is AE\ (r) + AEl (r) , different by about 
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100% from the London-only result for 3ao < r < 5ao. 
Because AEi(r) decreases with r as cxp(— 2r/ao), the 
London potential A-Ex(r) dominates at large r. But it 
is not until r = 5.78 a Q that \AE L (r)\ > 2|AJ5i(r)|. By 
then \AE L {r)\ < kT/5. 

What about second order beyond AEz,{r)l The cor- 
rection to the ground-state energy due to the perturba- 
tion W is 
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in which the sum is over all cigenstates \n) of the hamil- 
tonian Hq except |0). If the first state to contribute to 
the sum has energy Ei, then < E\ — Eq < E n — Eq 
or —1/(E\ — Eq) < —l/(E n — Eq). So the second-order 
change in the energy is bounded by 
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\(0\W\n}\ 2 < AE 2 (r) < 0. (4) 



The first states to contribute to the sum are those in 
which one of the electrons is in the n = 2 level, so E± = 
— (l/2)a 2 rac 2 (l + 1/4). But because the hydrogen wave 
functions fall off as exp(— rj(n do)), at long distances the 
first states to contribute are those with both electrons 
in the n = 2 level; they have Ei = -(l/2)a 2 mc 2 (l/4 + 
1/4). Without the energy denominators, the sum over 
intermediate states follows from completeness: 
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where (0|W|0) is the first-order result AEi(r), Eq.([2|). 

Evaluating the mean value of W 2 involves the integral 
G(x) = (2/tt) Jd 3 y(l/y 2 ) exp(— 2\y — x\), a function only 
of the length x = \x\ = \r\/an. We find 

G(x) = - [(1 + 2x) e- 2x Ei(2x) + (1 - 2x) e 2x Ei(2x)] 

(6) 

in which Ei and Ei are exponential integrals. G{x) con- 
tributes to the integral J(x) 

1 f°° r 
J{x) = — \ yG(y) (l + 2\x-y\) e' 2 ^ 
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that represents the square of the interaction of the two 
electrons. The product of electron-electron and electron- 
proton energies is described by the function K(x) 



K(x) = <Mln2-3 
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FIG. 1: The first-order van der Waals energy AEi(r), Eq.©, 
(solid, red) is nearly equal to the London result AEl(t-) (long 
dashes, green) when both are reliable and greater than kT 
(dotted, plum). The change in the sum AE(r) = AEi(r) + 
AE^r) (dashes, blue) as given by second-order perturbation 
theory, Eq.©, is 4 to 5 times greater than the London term 
for 3ao < r < 5an. The potential of Eq. (|15[) with a = 336 eV, 
b = 1.970Q 1 , c = 0.47«o 1 , d = 397ag eV, and e = 14093aJ 2 
(dashes, light blue) overlaps AE(r). 



in which 7 = 0.5772 ... is the Euler-Mascheroni con- 
stant. In terms of AEi, J, and K, the second-order 
correction AE 2 (r) to the energy at r = ao x becomes 
-(/a /e 2 )((0|^ 2 |0)-(0|^|0) 2 ) or 



AE 2 (r) = 




(9) 



+ — [G(x) + J(x)+K(x)] 



The fudge factor / = 1.083 = 6.499/6 is the one used by 
Pauling and Beach Q to estimate the effect of the energy 
denominators in the sum Eq.([3|). 

To second order in the perturbation W, the change 
in the energy of two hydrogen atoms at proton separa- 
tion r is the sum of the first-order term Eq.([2]) and the 
second-order term Eq.©, AE(r) = AEi(r) + AE 2 (r). 
Fig.CDshows AE(r) for 3 < r/a < 5 where AE is never 
even close to AE^. (To model non-reactive atoms, we 
have neglected the exchange terms that cause covalent 
bonding and only increase \AE(r)\.) 

What about additivity? Consider three hydrogen 
atoms at Ri, R 2 , and -R3. Again neglect exchange ef- 
fects and take the unperturbed ground state to be 

1,2,3) = |01,02,03) = \ls,Ri;ls,R 2 ;ls,R 3 ). (10) 
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Now the perturbation W is W 
which the pair potential Wij is 



W12 + W23 + W13 in 



\Rj-Ri\ \Rj + fj - Ri - n 



{Rj + fj-Ril \Rj-Rt 



(11) 



To first order in the perturbation W, the change in the en- 
ergy AEi = (1, 2, 3| W 12 + W 23 + W31 |1, 2, 3) is additive: 
AE 1 = <l,2|Wi a |l J 2) + (2,3|W23|2,3> + <l,3|Wi 3 |l,3). 

But at smaller separations the cross-terms do not van- 
ish. To second order in W, AE2 is not quite additive. 
The mean value of W12 W23, e.g., involves a non-zero 
sum over the intermediate states |l,n2,3) = |01,n2, 03), 

(1,2,3|VF 12 W 23 |1,2,3) = 

^(l,2 ) 3|Wi 2 |l,ri2 ) 3)(l,n2,3|W 2 3|l,2,3) 
n 

= £<l,2|W 12 |l,n2)(n2,3|W23|2,3>^0. (12) 

n 

Failure of additivity already occurs at second order. 

In the usual treatment, one expands the perturbation 
W in powers of 1/r keeping only terms as big as 1/r 6 
at large r and neglects terms like Eq. (fl~2|) . If R2 — Ri 
points in the 3-direction, then the leading term in the 
matrix element (1, 2| W12II, n2) involves [x\X2 + 2/12/2 — 
1z\Z2) l\R2 — -Ri| 3 which vanishes because of the spherical 
symmetry of the ground state of atom 1. This neglect of 
the first-order van der Waals force defers non-additivity 
to the third order of perturbation theory in the Axihod- 
Teller-Muto formulation @, This deferral, however, is 
probably inconsequential compared with the prior omis- 
sion of purely additive first-order interactions. 

Clearly there is room for imp rovement on the Lennard- 

jones @, i, EE, [HI, El EE El HE EE, EE El, EE EE m 

|22L |23| and harmonic forms 



V LJ (r) = \V(r s )\ 



V H (r) = V(r s ) + 
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(r-r s f d 2 V{r s ) 



2 dr 2 
What to do? We suggest using the form 

d 



V(r) 



- b r 
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) 



(13) 
(14) 

(15) 



The terms involving a, b, and c appear in Eq.([2]) for 
AEi(r), were proposed by Rydberg [24[ to incorporate 
spectroscopic data, but were largely ignored until re- 
cently [2E HE]. The constant d = Cq is the coefficient 
of the London tail. The new term e r -6 cures the Lon- 



don singularity; as r — > 0, d/(r 6 + er 



3 /e and 



V(r) — > a, finite. But as r — > oo, V(r) approaches the 
London term, V(r) — > —d/r 6 = —Cs/r 6 . 
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FIG. 2: The phenomenological potential (|15|l with a = 53.8 
eV, b = 2.99 A -1 , c = 2.453 A" 1 , d = C 6 = 3.884 eV A 6 , and 
e = 47.6 A 12 (solid, red) is finite, fits the RKR spectral points 
for molecular hydrogen (crosses, blue), and gives the correct 
London tail for r > 3 A. The harmonic potential Vh (dashes, 
green) is accurate only near its minimum. 



To test the ability of the potential V(r) to represent 
covalent bonds, we used Gnuplot [27J to choose a, b, c, d, 
and e so as to fit Eq. (|15[) to the experimental potentials 
of molecular IF , N 2 , and 2 obtained from spectroscopic 
data [lEIIEIaq] by the Rydberg 0-Klein (3l|- Rees^ 
(RKR) method. For H2, we set d equal to the London 
value, d = Ce = 6.5 mc 2 a 2 a®. Fig. §2fy shows that the 
potential of Eq. p^j) (solid, red) goes through the RKR 
points for H 2 (crosses, blue), while the harmonic poten- 
tial (|14p (dashed, green) is accurate only near its min- 
imum and is dreadful for r 3> r s . The fits of V(r) to 
the N 2 and 2 RKR points were equally good. The po- 
tential V(r) of Eq. (fT5)) represents covalent bonds better 
than the harmonic potential Vh of Eq. (fl~4|) (or Vlj). 

Can V(r) represent weak noncovalent bonds? Using 
Gnuplot, we fitted a, 6, c, and e in Eq JT5|) to calculated 
He2 [Si or empirical Ar 2 and Kr 2 |34| potentials, set- 
ting d equal to the respective London coefficient. Fig.([3|) 
shows that V(r) of Eq. p^|) (solid, red) fits the empirical 
points [34[ for Kr 2 (crosses, blue), but that the Lennard- 
Jones potential Vlj of Eq. (fT3|) (dashes, green) is too low 
for r > 4A. Fig. ((4|) shows that V(r) fits the empirical 
Kr 2 points, but that Vlj is much too hard for r < 4A. 
The fits for He2 and Ar 2 were equally good. Figs. ((2E1) 
would look even better with the London d treated as a 
free parameter. The potential V(r) of Eq. (fl~5|) represents 
weak noncovalent bonds better than Vlj or Vh- Vlj is 
unnaturally hard for r < r s ; it can wrongly reject Monte 
Carlo moves that put atoms slightly too close, frustrating 
M-C searches for the native states of macromolccules. 
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FIG. 3: The potential V of Kr-Kr attraction (Eg .1(15]) with 
a = 2296 eV, b = 2.5136 A" 1 , c = 0.2467 A -1 , d = 78.2146 
eVA 6 , and e = 407923 A 12 ) (solid, red) is finite, fits the 
Kr2 Aziz points (crosses, blue), and gives the correct London 
tail. When matched at the minimum, the Lennard- Jones form 
Vlj (Eq.([T3j with r s = 4.008 A and V(r s ) = -0.017338 eV) 
(dashes, green) is too low for r > 4 A. 
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FIG. 4: Positive potential V (Eq.|l5j) as in Fig. Q), Kr 2 Aziz 
points (crosses, blue), L-J form Vlj fEq. (|13|l as in Fig. <(3j ) . 

The van der Waals potential is often as large in first 
order as in second. The potential (|15p may be useful in 
macromolccular simulations that involve not only bonded 
but also non-bonded atoms. Given the differences be- 
tween Eq. (|15j) and the popular Lennard- Jones form, it 
would be worthwhile to examine the consequences of 
these differences in numerical simulations. 
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Saslow, and C. Schwieters for advice. P. J. Steinbach 
kindly hosted KC at NIH, where we used the Biowulf 
computers. 



* Electronic address: cahill@unm.edu 
' Electronic address: aparsegi@helix.nih.gov 
[1] L. D. Landau and E. M. Lifshitz, Quantum Mechanics 

(Pergamon, 1976), p. 341, 3rd ed. 
[2] B. M. Axilrod and E. Teller, J. Chem. Phys. 11, 299 
(1943). 

[3] Y. Muto, Phys.-Math. Soc. Jap. 17, 629 (1943). 

[4] R. Eisenschitz and F. London, Z. Physik 60, 491 (1930). 

[5] F. London, Z. Physik 63, 245 (1930). 

[6] L. Pauling and E. B. Wilson, Jr., Introduction to Quan- 
tum Mechanics (Dover, Mineola, NY, 1935), s. XIV-47a. 

[7] L. Pauling and J. Y. Beach, Phys. Rev. 47, 686 (1935). 

[8] H. Margenau, Rev. Mod. Phys. 11, 1 (1939). 

[9] R. Pastor, Ph.D. thesis, Harvard University (1984). 
[10] S. R. Atlas, Ph.D. thesis, Harvard University (1988). 
[11] T. A. Halgren, J. Am. Chem. Soc. 114, 7827 (1992). 
[12] W. J. Meath and D. J. Margoliash, J. Chem. Phys. 68, 
1426 (1978). 

[13] A. D. Buckingham, P. W. Fowler, and J. M. Hutson, 

Chem. Rev. 88, 963 (1988). 
[14] D. Ben-Amotz and D. R. Herschbach, J. Phys. Chem. 

94, 1038 (1990). 
[15] D. Ben-Amotz and K. G. Willis, J. Phys. Chem. 97, 7736 

(1993). 

[16] G. Chalasihski and M. Szcz§sniak, Malgorzata, Chem. 

Rev. 94, 1723 (1994). 
[17] B. Jeziorski, R. Moszynski, and K. Szalewicz, Chem. Rev. 

94, 1887 (1994). 
[18] F. B. van Duijneveldt, J. G. C. M. van Duijneveldt-van de 

Rijdt, and J. H. van Lenthe, Chem. Rev. 94, 1873 (1994). 
[19] D. Yin and A. D. MacKerell, Jr., J. Comp. Chem. 19, 

334 (1998). 

[20] T. C. Choy, Phys. Rev. A 62, 012506 (2000). 

[21] I. J. Chen, D. Yin, and A. D. MacKerell, Jr., J. Comp. 

Chem. 23, 199 (2001). 
[22] D. Ben-Amotz and K. G. Willis, J. Chem. Phys. 117, 

6590 (2002). 

[23] A. I. Volokitin and B. N. J. Persson, Phys. Rev. B 65, 

115419 (2002). 
[24] R. Rydberg, Z. Phys. 73, 376 (1931). 
[25] J. Ferrante, H. Schlosser, and J. R. Smith, Phys. Rev. A 

43, 3487 (1991). 
[26] J. N. Murrell, S. Carter, S. C. Farantos, P. Huxley, and 

A. J. C. Varandas, Molecular Potential Energy Functions 

(John Wiley & Sons, 1984), p. 9. 
[27] http://www.gnuplot.info/. 

[28] S. Weissman, J. T. Vanderslice, and R. Battino, J. Chem. 

Phys. 39, 2226 (1963). 
[29] A. Lofthus and P. H. Krupenie, J. Phys. Chem. Ref. Data 

6, 113 (1977). 

[30] P. H. Krupenie, J. Phys. Chem. Ref. Data 1, 423 (1972). 
[31] O. Klein, Z. Phys. 76, 226 (1932). 

[32] A. L. G. Rees, Proc. Phys. Soc. (London) 59, 998 (1947). 



5 



[33] J. B. Anderson, J. Chem. Phys. 115, 4546 (2001). 

[34] R. A. Aziz and M. J. Slaman, Mol. Phys. 58, 679 (1986). 



